
#include "common.h"

/*
 * CPU tranformation version
 *
 */
void corrcoefmatrix_tile(float *out, float *in, int n, int m, int blockIdx_x, int blockIdx_y, int threadIdx_x, int threadIdx_y) {

  float inAS[TILESIZE][TILESIZE];
  float inBS[TILESIZE][TILESIZE];
  int i,j,k;
  int tileBegin;
  float a_x,b_x;
  float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;
int tx,ty;

  i = blockIdx_y*TILESIZE + threadIdx_y;
  j = blockIdx_x*TILESIZE + threadIdx_x;

  sum_a = sum_a2 = sum_b = sum_b2 = sum_ab = 0;
  for (tileBegin=0; tileBegin < m; tileBegin += TILESIZE) { 

     // Load data into shared memory
     for (ty=0; ty < TILESIZE; ty++) {
        for (tx=0; tx < TILESIZE; tx++) {
           inAS[ty][tx] = in[(blockIdx_y*TILESIZE + ty)*m+tileBegin+tx];
           inBS[ty][tx] = in[(blockIdx_x*TILESIZE + ty)*m+tileBegin+tx]; 
        }
     }

     if (i <= j) {

     // Do partial calculation on each tile
     for (k=0; k < TILESIZE; k++) {
        if (tileBegin + k < m) {  // we should test if it is as good as padding
           a_x = inAS[threadIdx_y][k];
           b_x = inBS[threadIdx_x][k]; 
           sum_a += a_x;
           sum_a2 += a_x*a_x;
           sum_b += b_x;
           sum_b2 += b_x*b_x;
           sum_ab += a_x*b_x;
        } 
     } 
     }

  }
  if (i <= j) {
  corrcoef = (m*sum_ab - sum_a*sum_b)/sqrtf((m*sum_a2-sum_a*sum_a)*(m*sum_b2-sum_b*sum_b));
  out[i*n+j] = corrcoef;
  out[j*n+i] = corrcoef;
  }
}

void corrcoefmatrix_debug(float *out, float *in, int n, int m) {
  int blockIdx_x, blockIdx_y, threadIdx_x, threadIdx_y;
  int i,j;

  for (blockIdx_y=0; blockIdx_y < (n/TILESIZE); blockIdx_y++) {
     for (blockIdx_x=0; blockIdx_x <  (n/TILESIZE); blockIdx_x++) {
        for (threadIdx_y=0; threadIdx_y < TILESIZE; threadIdx_y++) {
           for (threadIdx_x=0; threadIdx_x < TILESIZE; threadIdx_x++) {
              corrcoefmatrix_tile(out,in,n,m,blockIdx_x,blockIdx_y,threadIdx_x,threadIdx_y);
           
           }
        }
     }
  }
}

